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In this paper a novel tool BioDiVinE for parallel analysis of biological models is presented. The 
tool allows analysis of biological models specified in terms of a set of chemical reactions. Chemical 
reactions are transformed into a system of multi-affine differential equations. BioDiVinE employs 
techniques for finite discrete abstraction of the continuous state space. At that level, parallel analysis 
algorithms based on model checking are provided. In the paper, the key tool features are described 
and their application is demonstrated by means of a case study. 



1 Introduction 

The central interest of systems biology is investigation of the response of the organism to environmental 
events (extra-cellular or intra-cellular signals). Even in procaryotic organisms, a single environment 
event causes a response induced by the interaction of several interwoven modules with complex dynamic 
behaviour, acting on rapidly different time scales. In higher organisms, these modules form large and 
complex interaction networks. For instance, a human cell contains in the order of 10,000 substances 
which are involved in 15,000 different types of reactions. This gives rise to a giant interaction network 
with complex positive and negative regulatory feedback loops. 

In order to deal with the complexity of living systems, experimental methods have to be supplemented 
with mathematical modelling and computer- supported analysis. One of the most critical limitations in 
applying current approaches to modelling and analysis is their pure scalability. Large models require 
powerful computational methods, the hardware infrastructure is available (clusters, GRID, multi-core 
computers), but the parallel (distributed) algorithms for model analysis are still under development. 

In this paper we present the tool BioDiVinE for parallel analysis of biological models. BioDiVinE 
considers the model in terms of chemical equations. The main features of the tool are the following: 

• user interface for specification of models in terms of chemical equations 

• formal representation of the model by means of multi-affine ODEs 

• setting initial conditions and parameters of the kinetics 
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• setting parameters of the discrete abstraction 

• graphical simulation of the discretized state space 

• model checking analysis. 

As an abstraction method, BioDiVinE adapts the rectangular abstraction approach of multi-affine 
ODEs mathematically introduced in Q and algorithmically tackled in [T71 [3 by means of a parallel 
on-the-fly state space generator. The character of abstraction provided by this discretization technique is 
conservative with respect to the most dynamic properties that are of interest. However, there is the natural 
effect of adding false-positive behaviour, in particular, the abstracted state space includes trajectories 
which are not legal in the original continuous model. 

The structure of the paper is the following. Section [2] gives a brief overview of the underlying 
abstraction technique and the model checking approach. Section [3] presents in step-by-step fashion the 
key features of current version of BioDiVinE. Section]?] provides a case study of employing BioDiVinE 
for analysis of a biological pathway responsible for ammonium transport in bacteria Escherichia Coli. 

1.1 Related Work 

In our previous work we have dealt with parallel model checking analysis of piece-wise affine ODE 
models JT4l . The method allows fully qualitative analysis, since the piece- wise affine approximation of 
the state space does not require numerical enumeration of the equations. Therefore that approach, in 
contrast to the presented one, is primarily devoted for models with unknown kinetic parameters. The 
price for this feature is higher time complexity of the state space generation. In particular, time appears 
there more critical than space while causing the parallel algorithms not to scale well. 

In the current version of BioDiVinE all the kinetic parameters are required to be numerically speci- 
fied. In such a situation there is an alternative possibility to do LTL model checking directly on numerical 
simulations 1 18 , 8 1. However, in the case of unknown initial conditions there appears the need to provide 
large-scale parameter scans resulting in huge number of simulations. On the contrary, the analysis con- 
ducted with BioDiVinE can be naturally generalised to arbitrary intervals of initial conditions by means 
of rectangular abstraction. 

Rectangular abstraction of dynamic systems has been employed in ifTTIl for reachability analysis. The 
method has been supported by experiments performed on a sequential implementation in Matlab. The 
provided experiments have showed that for models with 10 variables the reachability analysis ran out 
of memory after 2 hours of computation. This gave us the motivation to employ parallel algorithms. 
Moreover, we generalise the analysis method to full LTL model checking. 

There is another work that employs rectangular abstraction for dynamic systems (6). The framework 
is suitable for deterministic modelling of genetic regulatory networks. The rectangular abstraction relies 
on replacing S-shaped regulatory functions with piece-wise linear ramp functions. The partitioning of 
the state space is determined by parameters of the ramp functions. In contrast to that work, we consider 
directly general multi-affine systems with arbitrarily defined abstraction partitions. 

2 Background 

2.1 Modelling Approach 

The most widely used approach to modelling a system of biochemical reactions is the continuous ap- 
proach of ordinary differential equations (ODEs). We consider a special class of ODE systems in the 
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form x = f(x) where x = (x\ , . . . , x n ) is a vector of variables and / = (/i , . . . , f n ) : W 1 —> W 1 is a vector of 
multiaffine functions. A multiaffine function is a polynomial in the variables x\ , . . . , x n where the degree of 
any variable is at most 1. Variables X[ represent continuous concentrations of species. Multiaffine ODEs 
can express reactions in which the stoichiometry coefficients of all reactants are at most 1. The system 
of ODEs can be constructed directly from the stoichiometric matrix of the biochemical system |[T2ll . 



A + B^C 



Figure 1 : Example of a multi-affine system 

Multi-affine system is achieved from the system of biochemical reactions by employing the law of 
mass action. In Figure [T] there is an example of a simple biochemical system represented mathematically 
as a multi-affine system. If all the reactions are of the first-order, in particular, the number of reactants 
in each reaction is at most one, then the system falls into a specific subclass of dynamic systems - the 
resulting ODE system is affine. An example of an affine system is given in Figure [2] 

f=k 1 -A-k 3 -B + k 2 -C 
c §=k 1 -A-k 2 -C + k 3 -B 

Figure 2: Example of an affine system 

2.2 Abstraction Procedure 

The rectangular abstraction method employed in BioDiVinE relies on results of de Jong, Gouze fl31l 
and Belta, Habets [7j. Each variable is assigned a set of specific (arbitrarily defined) points, so-called 
thresholds, expressing concentration levels of special interest. This set contains two specific thresholds 
- the maximal and the minimal concentration level. Existence of these two thresholds comes from the 
biophysical fact that in any living organism each biochemical species cannot unboundedly increase its 
concentration. The intermediate thresholds then define a partition of the (bounded) continuous state 
space. The individual regions of the partition are called rectangles. An example of a partition is given in 
Figure [3] 

The partition of the system gives us directly the finite discrete abstraction of the dynamic system. 
In particular, BioDiVinE implements a (discrete) state space generator that constructs a finite automa- 
ton representing the rectangular abstraction of the system dynamics. Since the states of the automaton 
are made directly by the rectangles, the automaton is called rectangular abstraction transition system 
(RATS). The algorithm for the state space generator of RATS has been presented in [2j. The idea be- 
hind this algorithm is based on the results §T5\ H. The main point is that for each rectangle the exit 
faces are determined. The intuition is depicted in Figure [4] There is a transition from a rectangle to its 
neighbouring rectangle only if in the vector field considered in the shared face there is at least one vector 
whose particular component agrees with the direction of the transition. The important result is that in a 
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f = -k r A+k 2 -B \ \ j 

Ul g , 

f=h.A-k 2 -B - j 

thresholds on A: 0,5,6, 10 

thresholds on B: 0, 2, 3, 5 J LJ I 

A 

Figure 3: Example of a rectangular partition of a two-dimensional system 

















~ 1 





Figure 4: Intuition behind the construction of the rectangular abstraction transition system 

multi-affine system it suffices to consider only the vector field in the vertices of the face. In Figure |4j the 
exit faces of the central rectangle are emphasised by bold lines. In Figure [5] there is depicted the rect- 
angular abstraction transition system constructed for the affine system from Figure |3] It is known that 
the rectangular abstraction is an overapproximation with respect to trajectories of the original dynamic 
system. 

There is one specific issue when considering the time progress of the abstracted trajectories. If 
there exists a point in a rectangle from which there is no trajectory diverging out through some exit 
face, then there is a self-transition defined for the rectangle. In particular, this situation signifies an 
equilibrium inside the rectangle. Such a rectangle is called non-transient. For affine systems there is 
known a sufficient and necessary condition that characterises non-transient rectangles by the vector field 
in the vertices of the rectangle. However, for the case of multi-affine system, only the necessary condition 
is known. Hence, for multi-affine systems BioDiVinE can treat as non-transient some states which are 
not necessarily non-transient. 

2.3 Model Checking 

In the field of formal verification of software and hardware systems, model checking refers to the problem 
of automatically testing whether a simplified model of a system (a finite state automaton) meets a given 
specification. Specification is stated by means of a temporal logic formulae. In the setting of RATS, 
model checking can be used in two basic ways: 

1. to automatically detect presence of particular dynamics phenomena in the system 

2. to verify correctness of the model (i.e., checking whether some undesired property is exactly 
avoided) 

In the case of dynamic systems we use Linear Temporal Logic (LTL) (see JTOl for details on LTL 
syntax and semantics interpreted on automata). LTL can be also directly interpreted on trajectories 
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-si X- 



CD -CD 



CD 



CD -CD 



9= ^3=^3 



Figure 5: Example of a rectangular abstraction transition system 



CD— —CD -CD 



CC- -CV-CD 



Figure 6: A counterexample contradicting the property of reaching B > 3 in finite time 



of dynamic systems (see e.g. |[T3l for definition of the semantics). Given a dynamic system S with a 
particular initial state we can then say that S satisfies a formula <p, written S \= <p, only if the trajectory 
starting at the initial state satisfies (p. In the context of automata, LTL logic is interpreted universally 
provided that a formula <p is satisfied by the automaton A, written A\= (p, only if each execution of the 
automaton starting from any initial state satisfies (p. The following theorem characterises the relation 
between validity of <p in the rectangular abstraction automaton and in the original dynamic system. 



Theorem 2.1 Consider a dynamic system S and the associated RATS A. If A \= (p then S \= (p. 



The theorem states that when model checking of a particular property on a RATS returns true, we are 
sure that the property is satisfied in the original dynamic system. However, when the result is negative, 
the counterexample returned does not necessarily reflect any trajectory in the original system. 

The system in Figure [5] satisfies a formula FG(B < 3) expressing the temporal property stating that 
despite the choice of the initial state the system eventually stabilises at states where concentration of B 
is kept below 3. Now let us consider a formula F(B > 3) expressing the property that despite the initial 
settings, the concentration of B will eventually exceed the concentration level 3. In this case the model 
checking returns one of the counterexamples as emphasised in Figure [6] stating that if initially A < 5 and 
B < 3 then B is not increased while staying indefinitely long in the shaded state. 
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Figure 7: BioDiVinE Toolset Architecture 
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Figure 8: A biochemical model specified in BioDiVinE GUI 

3 BioDiVinE Tool 

BioDiVinE employs aggregate power of network-interconnected workstations (nodes) to analyse large- 
scale state transition systems whose exploration is beyond capabilities of sequential tools. System prop- 
erties can be specified either directly in Linear Temporal Logic (LTL) or alternatively as processes de- 
scribing undesired behaviour of systems under consideration (negative claim automata). From the algo- 
rithmic point of view, the tool implements a variety of novel parallel algorithms JUHl for cycle detection 
(LTL model checking). By these algorithms, the entire state space is uniformly split into partitions and 
every partition is distributed to a particular computing node. Each node is responsible for generating the 
respective state-space partition on-the-fly while storing visited states into the local memory. 

The state space generator constructs the rectangular abstraction transition system for a given multi- 
affine system. The scheme of the tool architecture is provided in Figure [7] Library-level components are 
responsible for constructing, managing and distributing the state space. They form the core of the tool. 
The tool provides two graphical user interface components SpecAff — allowing editing of biological 
models in terms of chemical reactions, and SimAff — allowing visualisation of the simulation results. 

The input (biochemical) model is given by the following data: 

• list of chemical species, 

• list of partitioning thresholds given for each species, 

• list of chemical reactions. 
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VARS:A,B,C 

EQ:dA = (-0.1)*A 

EQ:dB = 0.1*A + (-1)*B + 1*C 

EQ:dC = 0.1*A + 1*B + (-1)*C 



TRES : A : 0, 4, 6, 10 
TRES : B : 0, 2, 10 
TRES:C: 0, 2, 4, 10 



INIT: 4:6, 0:2, 2:4 



Figure 9: A multi-affine ODE model and its state space generated by BioDiVinE 
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Figure 10: Visualisation of the state space in the BioDiVinE GUI, states are projected onto the BC 
concentration plane 

The biochemical model together with the tested property and initial conditions are then automatically 
translated into a multi-affine system of ODEs forming the mathematical model that can be analysed by 
BioDiVinE algorithms. The mathematical model consists of the following data: 

• list of variables, 

• list of (multi-affine) ODEs, 

• list of partitioning thresholds given for each species, 

• list of initial rectangular subspaces (the union of these subspaces forms the initial condition), 

• Biichi automaton representing an LTL property (this data is not needed for simulation). 

An example of a simple three-species model representing three biochemical reactions A — > B + C, 
B ^ C is showed in Figure [I] The first reaction is performed at rate 0.1 s~ l . The second two reversible 
reactions are at rate 1 j" 1 . The respective mathematical model is showed in Figure [9] on the left in 
the textual .bio format. For each variable there is specified the equation as well as the list of real 
values representing individual threshold positions. The initial condition is defined in this particular case 
by a single rectangular subspace: A e (4, 6), 5 e (0,2), and C E (2,4). The state space generated for 
this setting is depicted in Figure [9] on the right. Figure [TO] demonstrates visualisation features of the 
BioDiVinE GUI. 

For model checking analysis, BioDiVinE relies on the parallel LTL model checking algorithms of 
the underlying DiVinE library [4j. A given LTL formula is translated into a Biichi automaton which 
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Figure 11: A visualisation of the state-spaces in BC projection with uniform (left) and automatic (right) 
gradual threshold refinement, manually augmented by numeric simulation trace from COPASI AH 



VARS:A,B,C 

EQ:dA = (-0.1)*A 
EQ:dB = 0.1*A + (-1)*B 
EQ:dC = 0.1*A + (1)*B + 



- 1*C 
(-1)*C 



TRES : A 
TRES : B 
TRES : C 



0, 4, 
0, 2, 
0, 2, 



6, 10 
10 

4, 10 



INIT: 4:6, 0:2, 2:4 



process LTL_propertyl { 

state ql, q2; 

init ql ; 

accept q2; 

trans 

ql -> q2 { guard B>4 kk B<4.5 

kk C>3 kk C<3.5; }, 

-> ql O, 
-> q2 {}; 



qi 
q2 

} 



system sync property LTL_propertyl ; 



process LTL_property2 { 

state ql, q2; 

init ql ; 

accept q2; 

trans 

ql -> q2 { 



guard B>4 kk 
kk C>5.5 



ql -> 
q2 -> 
} 



ql O, 
q2 O; 



B<4.5 
Ik C<6: 



system sync property LTL_property2 ; 



Figure 12: A multi-affine model extended with a never claim automaton for property 1 and property 2 



represents its negation. That way the automaton represents the never claim property. The automa- 
ton is automatically generated for an LTL formula and merged with the mathematical model by the 
divine . combine utility. An example of a model extended with a never claim property is showed in 
Figure [12] In particular, the first automaton LTL_propertyl specified in terms of the DiVinE language 
represents a never claim for the safety LTL formula G^(B >4A5<4.5AC>3AC<3.5) expressing 
that concentrations of species B and C will never enter the specified rectangular region. The unreacha- 
bility of a slightly different region is defined by the automaton for property 2. The results of the model- 
checking procedure are showed in Figure [13] Property 1 has been proved false by finding a run of the 
system visiting the specified region (states of run given as list), while property 2 has been proved true by 
extensively searching all the system runs and not finding any one that would cross the region specified in 
property 2. 

The choice of threshold values for each variable affects greatly the shape and size of the generated 
state-space. The refinement of a given partitioning — the addition of more thresholds to a set of former 
thresholds — may result in the unreachability of a part of a region reachable before. 

Since manual refinement of thresholds by adding numeric values can be tedious or unintuitive, two 
more advanced methods are available in BioDiVinE. The first method divides a given interval uniformly 
into subintervals of a given size, while the second method tries a more sophisticated automatic technique 
of dividing regions according to signs of the concentration derivatives inside these regions fTTll . The 



J. Barnat et.al. 



39 



Accepting cycle 



[pre-initial] 

[4(1),0(0),2(1)-PP:0] 

[4(1) ,2(1) ,2(1)-PP:0] 

[4(1),2(1),3(2)-PP:0] 

[4(1),4(2),3(2)-PP:1] 

[0(0),4(2),3(2)-PP:l] 

======= Cycle ======= 

[0(0),2(1),3(2)-PP:1] 
[0(0) ,2(1) ,2(1)-PP:1] 



No accepting cycle 



states: 33 

transitions: 81 

iterations: 1 

size of a state: 16 

size of appendix: 12 

cross transitions: 

all memory 56 . 5 MB 

time: 0.115177 s 



0: local states: 33 
0: local memory: 56.5 



Figure 13: Results of model-checking for property 1 (left), property 2 (middle), visualisation of reachable 
(green) and unreachable (red) regions 
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Figure 14: Scaling of model checking algorithms on a homogeneous cluster 



state-spaces resulting from the gradual application of both threshold refinement approaches are depicted 
in Figure [TT] It is important to mention that the overall size of the state-space depends exponentially on 
the number of thresholds for all species. However, in some cases the actual reachable subspace of the 
whole state-space may be only polynomial in the number of thresholds. 

For any multi-affine model extended with a never claim automaton as showed in Figure 12 the 
parallel model checking algorithms can be directly called. 

We have performed several experiments Q in order to show scaling of the algorithms when dis- 
tributed on several cluster nodes. Figures [T4| and [T5| show scaling of model checking conducted on a 
simple model of a reaction network representing a catalytic reaction scaled for different numbers of 
intermediate products. 



4 Case Study: Ammonium Transport in E. Coli 

In this section we present a case study conducted using the current version of BioDiVinE. Since the 
rectangular abstraction method of multi-affine systems implemented in BioDiVinE is a relatively new 
result of applied control theory, its application is still in the stage of experimentation. In fact, we are 
still unaware of any case studies that apply this method to real biological models. In this case study 
we focus on demonstrating the usability of rectangular abstraction to analysis of biological models. 
To this end, we consider a simple biological model that specifies ammonium transport from external 
environment into the cells of Escherichia Coli. This simplified model is based on a published model 
of the E. Coli ammonium assimilation system ||T9| which was developed under the EC-MOAN project 
(http : //www. ec-moan. org). The metabolic reactions and regulatory reactions in the original model 
were removed. 
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Figure 15: Scaling of model checking algorithms on a homogeneous cluster, X-axis shows number of 
cluster nodes, Y-axis gives speedup relative to system size and number of nodes 

4.1 Model Description 

E. coli can express membrane bound transport proteins for the transportion of small molecules from 
the environment into the cytoplasm at certain conditions. At normal ammonium concentration, the free 
diffusion of ammonia can provide enough flux for the growth requirement of nitrogen. When ammo- 
nium concentration is very low, E. coli cells express AmtB (an ammonia transporter) to complement the 
deficient diffusion process. Three molecules of AmtB (trimer) form a channel for the transportation of 
ammonium/ammonia. Protein structure analysis revealed that AmtB binds NHA + at the entrance gate of 
the channel, deprotonates it and conducts NH^ into the cytoplasm as illustrated in Figure [T6| (left) JT6l . 



At the periplasmic side of the channel there is a wider vestibule site capable of recruiting NH^ cations. 
The recruited cations are passed through the hydrophobic channel where the pKa of NH^ was shifted 
from 9.25 to below 6, thereby shifting the equilibrium toward the production of NH3. NH^ is finally 
released at the cytoplasmic gate and converted to NH^ because intracellular pH (7.5) is far below the 
pKa of NH+. 

In addition to the above mentioned AmtB mediated transport, the bidirectional free diffusion of the 
uncharged ammonia through the membrane is also included in the simplified model. The intracellular 
NH^ is then metabolised by Glutamine Synthetase (GS). The whole model is depicted in Figure 



16 



(right). The external ammonium is represented in the uncharged and charged forms denoted NH^ex and 
NH^ex. Analogously, the internal ammonium forms are denoted NH^in and NH^in. The biochemical 
model that combines AmtB transport with NH^ diffusion is given in Table [T] The kinetic parameters are 
based on the values in the original model. 

By employing the law of mass action kinetics the reaction network is transformed into the set of 
multi-affine ODEs as listed in Table [2j Since we are especially interested in how the concentrations 
of internal ammonium change with respect to the external ammonium concentrations, we employ the 
following simplifications: 

• We do not consider the dynamics of the external ammonium forms, thus we take NH^ex and 
NH^ex as constants (the input parameters for the analysis). 
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Figure 16: E. Coli ammonium transport mechanism and the respective pathway 



AmtB + NH 4 ex ^-h AmtB : NH 4 
AmtB : NH 4 ^ AmtB : NH 3 + H ex 

AmtB : NH 3 ^ AmtB + NH 3 in 

k 5 



NH 3 in + H l 



NH 4 in 

k 6 k 7 



NH 4 in 



NH 3 ex &-^> NH 3 in 



Jfei =5-10 8 ,£ 2 = 5-10 3 
yt 3 = 50 

>fc 4 = 50 
yfe 5 = 80 

^ 6 = 1 . 10 15 ,yfc 7 = 5.62- 10 5 
yfc 8 = £ 9 = 1.4-10 4 



Table 1 : Biochemical model of ammonium transport 



• We assume constant intracellular pH (7.5) and extracellular pH (7.0), thus H ex and H in are calcu- 
lated to be 3 • 10~ 8 and 10~ 7 . Based on the extracellular pH and the total ammonium/ammonia 
concentration, concentrations of NH 3 ex and NH^ex can be calculated. 

Without the loss of correctness, we simplify the notation of the cation NH^ as NH 4 . 



4.2 Model Analysis 

From the essence of biophysical laws, it is clear that the maximal reachable concentration level accu- 
mulated in the internal ammonium forms directly depends on the ammonium sources available in the 
environment. However, it is not directly clear what particular maximal level of internal ammonium is 
achievable at given amount of external ammonium (distributed into the two forms). In the analysis we 
have focused on just this phenomenon. More precisely, the problem to solve has been to analyse how 



= -k x • AmtB • NH 4 ex + k 2 • AmtB : NH 4 + k 4 • AmtB : NH 3 

dAmtB:NH 3 = ^ . ^ _ ^ AmB . ^ 

dAmtB:NH 4 = ^ . . jjjj^ _ ^ . AmB . jjjj^ _ ^ . AmB . N jj a 

^^=k 4 -AmtB:NH 3 -k 7 -NH 3 in + k 6 -NH 4 in 
= k 5 'NH 4 in + k 7 NH 3 in H in -k 6 NH 4 in 



Table 2: Mathematical model of ammonium transport 
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Partition 


AmtB AmtB\NH 3 AmtB : NH 4 NH 3 in NH 4 in 


(1) 

(2) 


7 3 3 6 4 
7 9 3 8 26 



Table 3: Numbers of thresholds in partitions (1) and (2) 



the setting of the model parameters NH 3 ex and NH^ex affects the maximal concentration level of NH 3 in 
and NH^in reachable from given initial conditions. 

During discussions with biologists we have found out that there is currently not available any in 
vitro measurement of the Am^-concentration (and also the concentrations of dimers AmtB : NH 3 and 
AmtB : NH4). Hence there has appeared the need to analyse the model with uncertain initial conditions. 
Such a setting fits the current features of BioDiVinE, especially the rectangular abstraction that naturally 
abstracts away the exact concentration values up to the intervals between certain concentration levels. 

4.2.1 Maximal reachable levels of internal ammonium forms 

At first, we have considered the (abstracted) initial condition to be set to the following intervals between 
concentration values: 

AmtB e (0, 1 • 1(T 5 ), AmtB : NH 3 e (0, 1 • 1(T 5 ), AmtB : NH 4 e (0, 1 • 1(T 5 ), 
NH 3 in e (1 • l" 6 , 1.1 • 1(T 6 ), NH 4 in e (2- 1(T 6 ,2.1 • 1(T 6 ) 

Note that the upper bounds as well as the initial intervals of internal ammonium forms have been set 
with respect to the available data obtained from the literature. 

We have considered two rectangular partitions. The partition (1) has been set basically according 
to the initial conditions. The partition (2) has been obtained by running one iteration of the automatic 
threshold refinement procedure to partition (1). Numbers of thresholds per each variable are given in 
Table El 

We have conducted several model checking experiments in order to determine the maximal reachable 
concentration levels of NH 3 in and NH^in. In particular, we have searched for the lowest a satisfying 
the property G(NH 3 in < a) and the lowest j8 satisfying G(NH 4 in < j8). The property Gp requires 
that all paths available in the rectangular abstraction from the states specified by the initial condition 
must satisfy the given proposition p at every state. Note that if the model checking method finds the 
property Gp false in the model, it also returns a counterexample for that. The counterexample satisfies 
the negation of the checked formula, which is in this case F^p. Interpreting this observation intuitively 
for the given formulae, we use model checking to find a path on which the species NH 3 in (resp. NH 4 in) 
exceeds the level a (resp. /3). 

We did not want to get precise values of a, j8, but we rather wanted to get their good approximation. 
At the starting point, we substituted for a (resp. j8) the upper initial bounds of the respective variables. 
Then we found the requested values by iteratively increasing and decreasing a (resp. j8). The obtained 
results are summarised in Table 01 

The results have shown that NH 3 in does not exceed its initial level no matter how the external am- 
monium is distributed between NH 3 ex and NH^ex. The upper bound concentration considered for both 
NH 3 ex and NHA + ex has been set to 1 • 10~ 5 which goes with the typical level of concentration of the gas 
in the environment. 

In the case of NH 4 in we have found that the upper bound to maximal reachable level is in the in- 
terval /3 E (5.3 • 10 _4 ,5.4 • 10~ 4 ). Since the counterexample achieved can be a spurious one due to the 
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Partition (1) 


Partition (2) 


a 


GyNH^in < a) 


# states 


Time (# nodes) 


v*(NH'$in < OC) 


# states 


lime (# nodes) 


1 . 1 • 1U 


true 


1 CiQ 1 
lUol 


U.Jo s (1 J 


true 


1 < 1 n5 
l.j • 1U 


l.y S (lo) 


P 


G(NH 4 in<l3) 


# states 


Time (# nodes) 


G(NH 4 in<l3) 


# states 


Time (# nodes) 


1 • 10" 3 


true 


2161 


0.45 s (1) 


true 


1.6- 10 5 


2 s (18) 


5 • 10" 4 


false 


4753 


1.9 s (1) 


false 


2.7 -10 5 


3 s (18) 


6 10" 4 


true 


2161 


0.43 s (1) 


true 


1.5 -10 5 


1.8 s (18) 


5.4 -lO" 4 


true 


1441 


0.27 s (1) 


true 


2.1 -10 5 


4.2 s (18) 


5.3 -lO" 4 


false 


3421 


1.2 s (1) 


false 


2.7 -10 5 


2.2 s (18) 



Table 4: Experiments on detecting maximal reachable levels of internal ammonium forms 



Partition (1) 


Partition (2) 


NH 3 ex 


9 


# states 


Time (# nodes) 




# states 


Time (# nodes) 


19.5- 10- 4 

19.6 - 10- 4 


true 
false 


901 
1261 


0.22 s(l) 
0.6 s(l) 


true 
false 


1.4 -10 5 
3.4- 10 5 


1.9 s (36) 
5.9 s (36) 



Table 5: Experiments on detecting NH^ex levels affecting maximal reachable NH^in 



overapproximating abstraction, the exact maximal reachable value may be lower. To this end we have 
conducted several numerical simulations which give us the argument that our estimation of j8 is plausible. 

4.2.2 Dependence of stable internal ammonium on changes in external conditions 

In the second experiment, we have focused on determining how much external ammonium have to be 
increased in particular form in order to stimulate NH3in to exceed the considered initial level. The setting 
of partitions and initial conditions has been considered the same as in the previous experiments. 

First, we have varied the constant amount of NH^ex to find at which level of NH^ex the maximal 
reachable level of NH^in is affected. More precisely, we have observed for which setting of NH^ex 
the property (p = G{NH^in < 1.1 • 10 -6 ) is true and for which it is not. The relevant experiments are 
summarised in Table [5] We have found out that if NH^ex is set to 19.6 • 10~ 4 or higher level then 
NH^in increases above the upper initial bound 1.1 • 10~ 6 . The counterexamples returned again agree 
with numerical simulations. 

Finally, we have varied the amount of NH^ex in order to find at which level of NH^ex the maximal 
reachable level of NH^in is affected. The results presented in Table [6] give us the conclusion that despite 
the level of NH^ex (checked up to 10 12 ), the maximal level reached by NH^in remains the same. In 
particular, NH^in does not exceed the initial bounds. 



Partition (1) 


NH 4 ex 


<P 


# states 


Time (# nodes) 


1 

1 - 10 12 


true 
true 


901 
901 


0.25 s (1) 
0.25 s (1) 



Table 6: Experiments on detecting NH^ex levels affecting maximal reachable NH^in 
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4.2.3 Performance 

All the experiments have been performed on a homogeneous cluster allowing computation on up-to 22 
nodes each equipped with 16GB of RAM and a quad-core processor Intel Xeon 2GHz. The model we 
have dealt with contained 5 dynamic variables and 2 constants. With partition (1) the generated state 
space had maximally 4753 states and with partition (2) maximally 3.4 • 10 5 . When trying to run one 
more step of automatic partition refinement, the number of thresholds exceeded the memory reserved for 
storing of the mathematical model. Note that the model has to be stored in the memory of each node. 
Only the state space is distributed over the cluster. 

5 Conclusion 

In this paper, we have presented the current version of BioDiVinE tool which implements rectangular 
abstraction of continuous models of biochemical reaction systems. The tool provides the framework for 
specification and analysis of biochemical systems. The supported analysis technique is based on the 
model checking method. Linear temporal logic is used for encoding of the properties to be observed on 
abstracted systems. 

The tool provides parallel model checking algorithms that allow fast response times of the analysis. 
We have provided a case study on which the key features of the tool are demonstrated. The case study 
has showed that the tool can be used for quickly getting the approximation of maximal reachable con- 
centration levels of individual species in the model. In general, we have analysed the model with respect 
to the set of safety properties which are technically tackled by construction of the state space reachable 
from the given initial states. We have found out that the main advantage of the rectangular abstraction is 
the possibility to analyse the system with uncertain initial conditions. 

The current drawback of the abstraction method is strong overapproximation of non-transient states 
in multi-affine systems. In consequence, analysis of liveness kind of properties (e.g., oscillations, insta- 
bility) is infeasible because of large amount of spurious counterexamples that come from false identifica- 
tion of non-transient states. However, this is not the case for affine systems on which liveness properties 
can be checked without these problems. Since the analysed systems are typically non-affine, we can still 
employ the liveness checking on their linearised approximations. However, by the linearisation process 
the precision of the analysis is significantly reduced. Improving the tool in these aspects requires further 
research. 

In general, we leave for future work the development of methods for identification of spurious paths. 
We think that one of the promising directions in using the discrete abstractions for analysis of biologi- 
cal models is employing the model-checking-based analysis for extensive exploration of properties. In 
particular, instead of returning only one path, the model checker should provide a set of paths. In this 
directions we aim to continue the research. 
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